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A SPACE TIME VARYING COEFFICIENT MODEL: THE EQUITY 
OF SERVICE ACCESSIBILITY^ 

By Nicoleta Serban 

Georgia Institute of Technology 

Research in examining the equity of service accessibility has emer- 
ged as economic and social equity advocates recognized that where 
people live influences their opportunities for economic development, 
access to quality health care and political participation. In this re- 
search paper service accessibility equity is concerned with where and 
when services have been and are accessed by different groups of peo- 
ple, identified by location or underlying socioeconomic variables. Us- 
ing new statistical methods for modeling spatial-temporal data, this 
paper estimates demographic association patterns to financial service 
accessibility varying over a large geographic area (Georgia) and over 
a period of 13 years. The underlying model is a space-time varying 
coefficient model including both separable space and time varying co- 
efficients and space-time interaction terms. The model is extended to 
a multilevel response where the varying coefficients account for both 
the within- and between-variability. We introduce an inference pro- 
cedure for assessing the shape of the varying regression coefficients 
using confidence bands. 

1. Introduction. 

"A home is more than a shelter — when located in a community with resources 
and amenities it is a critical determinant of opportunity." [Blackwell and Fox 
(2004)] 

Service accessibility equity is the study of systematic disparities in a pop- 
ulation's access to services that are considered fundamental in fostering eco- 
nomic development, improving wellness and enhancing the general quality 
of life of a population within a given geographic area. Examples of such ser- 
vices are health care, education, healthy food, financial services and others. 
Accessibility is measured as utilization-scaled travel cost of a community U 
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to the nearby sites in a service network consisting of multiple service sites 
geographically distributed: S = {si, . . . , A common utilization measure 
is the population rate within the community and its surroundings [Marsh 
and Schilling (1994)]. In this paper the utilization is measured as the popu- 
lation rate divided by the service rate to account for the service availability 
for each population unit. 

One challenge in measuring service accessibility is defining the travel cost 
for the residents in a community to access the sites in the service network. 
In the research works so far, the travel cost is calculated as the average 
or minimum distance between the centroid of the region U and the nearby 
sites in the service network [Lovett et al. (2002); Talen (1997); Talen (2001)]. 
However, communities occupy uneven geographic areas varying in size, and, 
therefore, their simplified representation by their centroids is restrictive. In 
this research paper, we instead represent a community by a sample of loca- 
tions in the neighborhood U ,ui,. . . , ub S U ^ and compute the street-network 
distances from these sample locations to the service network. Furthermore, 
the travel cost at each sample location ui, is measured as a summary of the 
travel distances, {c?(?X6, Si)}i=i^...^„. 

Combining the two ideas discussed above, utilization- adjustment of the 
travel cost and representation of a community by a series of sample spatial 
points, we evaluate the accessibility of a neighborhood to a service network 
in year t using 



where C{ui„t) is the travel cost at the sample location Uj, measured as the 
average street-network distance to the closest Q service sites available at 
time t (in our study, (5 = 3), W{ub,t) is the utilization adjustment factor at 
location Ub and /3 is a distance utility parameter. We estimate /3 by robust 
linear regression: log{W{ub,t)) ~ —log{C{ub,t)). 

Dividing the geographic space into contiguous spatial units Us, s = l, . . . , S, 
where each spatial unit corresponds to a neighborhood (e.g., census tract), 
the accessibility measure (utilization-adjusted travel cost) varies across the 
geographic space and time; Y{Us,t) = Y{s,t) defines the space-time varying 
accessibility process. Moreover, there are multiple providers in the service 
network, the accessibility process has an intrinsic multilevel structure. Under 
this multilevel structure, let Yp{s,t) be the accessibility of the community Ug 
to the sites of the pth provider for p = 1, . . . ,P, where P is the number of 
service providers. 

This research paper focuses on measuring and estimating spatial-temporal 
patterns in the association between demographic variables (including race, 
ethnicity and income) and service accessibility. Specific questions that will be 
addressed within this study are as follows: Is service accessibility equitable 





A SPACE-TIME VARYING COEFFICIENT MODEL 



3 



across population groups varying in ethnicity and income? Do service distri- 
bution inequities vary across regions and time? Are there service providers 
that provide a more equitable distribution of their services than others? 
What is the most common demographic feature associated with inequities? 

To evaluate the equity of service accessibility with respect to various pop- 
ulation groups over a period of time and within a large geographic space, 
we propose to estimate the space-time varying association of the accessibil- 
ity measure jointly over a series of demographic variables. In this context, 
weak associations or the absence of systematic disparities in service access 
are interpreted as service accessibility equity. One challenge of this associa- 
tion analysis is simultaneous estimation of the association patterns since the 
goal is to assess both the equity with respect to various ethnicity and race 
demographic variables controlling for income and the equity with respect to 
income controlling for ethnicity and race. A second challenge is that services 
are delivered within a multilevel network - multiple providers which deliver 
across multiple service sites. 

Many existing studies have analyzed service accessibility for different 
groups of people identified by underlying socioeconomic variables, but they 
are limited to small geographic areas and to only one year of data [Graves 
(2003); Larson (2003); Poweh et al. (2007); Small and McDermott (2006); 
Talen (2001); Talen and Ansehn (1998); Zenk et al. (2005)]. Commonly 
employed statistical procedures include regression methods assuming inde- 
pendence between service sites. Exploratory studies rely on graphical diag- 
nostics but not on statistical inference, which can be used to make informed 
decisions. Although the methods applied to the existing studies have useful- 
ness for some research questions, a spatial-temporal multivariate analysis of 
data with a multilevel structure requires new statistical methods which are 
rigorous, take into account the dependence in the data, and implementable, 
apply to real data complexity. 

To this end, we introduce a space-time (multilevel) model which allows 
estimation of space-time varying association patterns of a set of functional 
predictors (e.g., demographic variables) to a functional response, in our 
case study, the accessibility process. The modeling procedure introduced 
in this paper falls under a more general framework: varying-coefficient mod- 
els. These models have been applied to longitudinal data to estimate time- 
dependent effects on a response variable [Assuncao (2003); Fan and Zhang 
(2000); Hastie and Tibshirani (1993); Hoover et al. (1998); Huang, Wu and 
Zhou (2002); Wu and Liang (2004); Zhang (2004)]. Waller et al. (2007) 
review existing models to explore space-varying regressions and propose 
a Bayesian procedure. Gelfand et al. (2003) briefly mention the extension 
of their proposed Bayesian space-varying model to separable space-time 
varying coefficient models with a warning on its computational challenges. 
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Space-time separability greatly simplifies the problem by reducing the com- 
putational effort; however, it is a restrictive assumption since it implies that 
dependence attenuates in a multiplicative/additive manner across space and 
time. Therefore, extension of the Bayesian varying coefficient model to more 
complex modeling (e.g., space-time interaction) requires expensive compu- 
tations which may be prohibitive for densely sampled spatial domains. 

Our methodological contribution is three- fold. First, we propose a space- 
time varying coefficient model that takes into account the interaction be- 
tween time and space in a computationally efficient manner. To overcome 
the computational complexity due to operations with a large dependence 
matrix, we use a low-rank approximation to the space-time coefficient pro- 
cesses using radial basis of functions [Ruppert, Wand and Carroll (2003)]; 
this approach enables estimation of the space-time varying coefficient model 
for densely observed space and/or time domains. 

Second, we extend this model to multilevel data, resulting in a multi- 
level varying coefficient model. A few recent works have considered the 
study of multilevel functional models [Baladandayuthapani et al. (2008); 
Crainiceanu, Staicu and Di (2009); Di et al. (2009); Morris and Carroll 
(2006); Morris et al. (2003); Rice and Wu (2001); Staicu, Crainiceanu and 
Carroll (2010); Wu and Zhang (2002)]. In the related research, models of 
multilevel functional data have been applied to functional responses where 
the predictor is a fixed variate, commonly time, and, more recently, they 
have been extended to functional predictors but scalar responses. In this pa- 
per, the multilevel functional model applies to both functional response and 
functional predictors and it extends to the more difficult setting when the 
functionality is with respect to space and time. Challenges in estimating such 
a complex model include nonidentifiability and computational efficiency. We 
overcome the identifiability problem by using a knots-based kernel decom- 
position with a different set of knots across the model coefficients. We use 
penalized splines for computational efficiency in adapting to the smoothness 
in the space-time varying coefficients [Ruppert, Wand and Carroll (2003)]. 

Third, we introduce an inference procedure to assess the shape of the 
space-time varying coefficients. Generally, a common approach for identify- 
ing the shape of a regression function is hypothesis testing. However, for our 
model, hypothesis testing will require multiple tests for deciding whether its 
shape is nonlinear, linear or constant as a function of space or/and time. In 
this paper, we discuss an inference procedure for assessing the shape of the 
varying regression coefficients using confidence bands. 

The rest of the paper is organized as follows. In Sections 2 and 3 we present 
the space-time varying coefficient model as well as its extension to multi- 
level data along with the estimation and inference procedures. In Section 4 
we present the application of the models introduced in this paper to evalu- 
ate the equity of financial service accessibility in Georgia. We first describe 
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the data resources followed by the discussion of our results and findings. 
Section 5 concludes the paper. Some technical details are deferred to the 
supplemental material [Serban (2011)], which also provides complemental 
graphical descriptions of our analysis of the equity of service accessibility. 

2. Space time varying coefficient model. 

2.1. The model. In this section we introduce a space-time varying co- 
efficient model for estimating the relationships between the accessibility 
process and a series of demographic variables varying in time and space. 
The observed data are {Yij , {Xr^ij , r = 1,...,R}), where Yij =Y{ti,Sj) is 
the response variable and Xj-^ij = Xr{ti,Sj) a set of covariates observed 
at location sj = (sji, Sj2), j = ^, ■ ■ ■ , S , and time tj, i = 1, . . . , T, such that 

E[yij|X] ='yi{ti,Sj)Xi^ij-\ \-^Riti,Sj)XR^ij where 7r(t,s) for r = 1, . . . , i? 

are smooth coefficient functions. Note that not all covariates need to vary 
in both time and space; the modeling procedure allows for various predictor 
forms (scalar, varying in time, varying in space or both). For example, in 
our model implementation we take = 1 and, therefore, 7i(,s,t) is the 
intercept coefficient. 

In this paper we decompose the regression coefficients into separable space 
and time global effects along with space-time deviations from the global 
effects which are intrinsically local and account for the interaction between 
space and time: 

Mr Nr 

Jr{t,s)=ar{t)+Pr{s) + X] ^^''""^t^'^pd* " '^m ^ I )^^p(ll « " II ) ' 

m=l n=l 

We decompose the global coefficient functions using the radial spline basis 
[Ruppert, Wand and Carroll (2003)], 

Mr 

(2.1) ar{t) = Trfl + Tr,lt+ ^ U^.m^^temp ( |t " '^m^D, 

771=1 

Nr 

(2.2) f3r{si,S2) = 6r,0 + 6r,llSi + 5r,12S2 + '^Vr,nKsp{\\s - K^f^||). 

n=l 

(T) 

In these decompositions Ktcmp{t) is a temporal kernel whereas Km , m = 
1, . . . , Mr, are knots covering the time domain, and Ksp{s) is a spatial kernel 

(s) 

whereas Kn , n = 1, . . . , A^^, are knots covering the space domain. 

Importantly, although the kernel of the space-time interaction coefficient 
is separable in time and space, the decomposition is not. One advantage of 
using this kernel decomposition is that it allows decomposition of the design 
matrix as a Kronecker product, which, in turn, will ease the computations in 
the estimation procedure. We derive the Kronecker product decomposition 
in the Supplemental Material 1 of this paper. 
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In the semiparametric literature a common kernel function is the radial 
spline kernel function defined for d-dimensional domains [Nychka and Saltz- 
man (1998)]. Bivariate smoothing based on radial basis functions has the 
advantage of being rotational invariant, which is important in geographical 
smoothing. For two-dimensional domains, the function Ksp{-) could be re- 
placed by any other covariance function [Cressie (1993)], for example, the 
Mater n covariance function. 

Knots for one-dimensional spaces are commonly set to the sample quan- 
tiles of the observation points, whereas knots for two-dimensional spaces are 
commonly selected using the space-filling algorithm [Nychka and Saltzman 
(1998)], which is based on minimax design, or /c-nearest neighbor clustering 
algorithms. In this paper we implement these standard methods to select 
the number of knots. 

2.2. Estimation. We choose a method for estimating the model described 
in the previous section from among several candidate procedures. One mod- 
eling approach is smoothing splines [Wahba (1990)], which assumes that the 
number of knots is equal to the number of observation design points {Mr = T 
and Nr = S) and controls the smoothness of the coefficient by penalizing the 
influence of the coefficients Ur^m, m= 1, . . . , M^, and Vr^n, n = 1, . . . ,Nr, us- 
ing a penalty function. One primary drawback of this estimation procedure 
is its computational aspect. A less computational approach is regression 
splines [Wahba (1990)], in which a small number of knots are used (M,. <C T 
and Nr <^ S). This reduces to selection of the optimal numbers of knots, 
which can be computationally expensive in the context of our model since 
it requires solving a multidimensional optimization problem. The smooth- 
ness levels of the regression functions differ from one covariate to another 
and, therefore, we need to optimally identify {Mr,Nr)s for r = 1, . . . ,R. In 
addition, this approach introduces modeling bias. 

An alternative approach to optimal knots selection is to assume equal 
number of spatial knots [Nr = N for r = 1, . . . ,R) and equal number of tem- 
poral knots {Mr = M for r = 1, . . . , i?) with N and M sufficiently large such 
that the modeling bias is small [Li and Ruppert (2008)], but, similarly to 
smoothing splines, impose constraints on the coefficients Ur^m-, 'Vr,n, and fr,nm 
as follows: 

M N M N 

E<n.<CP, T.<n<Ci'\ EE<n.n<Cr 

m=l n=l m=ln=l 

or, equivalently, estimate the coefficients using penalized regression 

\\h(yk,ij) -ll{ti,Sj)Xi^ij -fr{ti,Sj)Xr,ij\\'^ 

R 

+ ^{A(^)u,.u; + A(^V,.v; + Xrl^rK}^ 
r=l 
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where 

Ur = {Ur,l,...,Ur,N), = {Vr,l, . . . ,Vr,M) and 

l/f. = {Vr,nm}n=l,...,N,m=l,...,M- 
(T) (S) 

Moreover, the parameters Ar , Ar and A^ are penalties controlhng the 
smoothness level of the regression coefficients. This approach is often re- 
ferred to as penalized splines [Ruppert, Wand and Carroll (2003)]. Conse- 
quently, selection of the number of knots reduces to selection of the penalty 
parameters, which, in turn, is a multidimensional optimization problem. 

In the semiparametric regression literature the problem of selecting the 
penalties, and implicitly of the the number of knots, is overcome by solving 
an equivalent mixed effects regression problem where u^, and are 
random effects, specifically, ~ N{0, {ar^^)'^In)-, ~ A^(0, {ar^^)'^Im) and 
Ur ~ N[0,a'^lnm)- Under the mixed effects model, the penalties are 

AI. ^ — TT^^T , a1 ^ — 7777 and A7- — ^ . 

We therefore estimate the model parameters using a mixed effects model to 
circumvent the difficulty of selecting the penalty parameters, or, implicitly, 
the number of knots. 

Based on the mixed-effects model formulation, denote the vector of the 
fixed effects 

= [n,o n,i (5i,ii (5i,i2 ••• tr^ tr^i Sr^ Sr^h Sr^i2] 

with identifiability constraints 5rfi = for r = 1, . . . , R. The vector of random 
effects is 

U=[ui Vi Vi ■■■ UR VR Vr]. 

The corresponding design matrices X and Z are 

X = [Xi{ti,Sj) ■■■ X/j(ti,Sj)]j^^ 7,^^;^ 5 

with Xr{ti, Sj) = Xriti,Sj)[l ti 1 Sij S2j] , 
Z=[Zi{ti,Sj) ■■■ 'ZR{ti,Sj)]i^-^_^j,j^-i^_g 
with Zr{U,Sj)=Xr{k,Sj) [Ktcmpdi* " «^^^|) Ksp{\\Si - K^f^W) 

-^tcmp ( i l^rn ^ I ) Kgp ( I [ S j Kn | [ ) ] • 

The model in the matrix form becomes E[y |^] = X& + ZXJ, which is equiv- 
alent to a linear mixed model. 
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2.3. Inference. In this section we discuss alternative methods for making 
inference on the shape of the regression coefficients entering the space-time 
varying coefficient model described in Section 2.1. Specifically, we discuss 
a procedure for evaluating the shape (constant vs. linear vs. nonlinear) of 
the temporal and spatial global coefficients and a procedure for testing the 
significance of the space-time interaction. 

Shape evaluation. In this section we discuss a novel procedure for shape 
evaluation of the temporal regression coefficients. A similar procedure ap- 
plies to spatial regression coefficients. Specifically, each temporal regression 
coefficient can take various shapes, for example, constant [a{t) = tq], linear 
[a{t) = To + Tit] or nonlinear. 

In varying-coefficient models, the common procedure for assessing the 
shape of the coefficients is hypothesis testing. For example, the hypothesis 
test for linearity is equivalent to Hq : ct^ = vs. Ha :a^> where o"^ is the 
variance of the random effects Um under the mixed effects model. The com- 
mon approach is a likelihood ratio testing (LRT) procedure. Crainiceanu et 
al. (2005) developed a LRT by taking advantage of the existing research in 
hypothesis testing for zero variance in linear mixed-effects (LME) models. 
Liang, Wu and Carroll (2003) tested for linearity of nonparametric functions 
using a Cramer-von Mises statistic. 

Although there are several competitive approaches for testing for linearity 
of the regression coefficients, because we need to test sequential hypotheses 
to decide about the shape of a coefficient and because we often have a large 
number of predictors that enter the space-time varying coefficient model, we 
instead propose identifying the shape of the coefficients using simultaneous 
confidence bands. If CB^ is a 1 — 7 confidence band for the coefficient a{t), 
then P{a{t) € CB^, t € T) > 1 — 7 where T is the time domain. The deriva- 
tion of the joint confidence bands is presented in the Supplemental Material 2 
of this paper. 

Many authors have noted that using confidence intervals has a series of 
advantages over the conventional hypothesis testing [Sim and Reid (1999)]. 
Confidence intervals cannot only be used to test a hypothesis, but also to pro- 
vide additional information on the variability of an observed sample statistic 
and on its probable relationship to the value of this statistic in the popula- 
tion from which the sample was drawn. 

Figure 1 depicts examples of three different one-dimensional shapes along 
with their confidence bands (CB). We define "constant" shape if there exists 
a constant line that falls within the confidence bands. Similarly, we define 
"linear" shape if there exists a linear function that falls within the confidence 
bands. When searching for a line L within the confidence bands, it suffices 
to search for linear functions between the convex hull of the upper level of 
the confidence band and the convex hull of the lower level of the confidence 
band. Although the coverage of the shape test hypothesis is maintained when 
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Fig. 1. Examples of nonlinear, linear, constant effects. 

using confidence bands, the power may be reduced, as we point out in the 
simulation study included in Supplemental Material 4. 

According to the result below, accepting the null hypothesis at a signifi- 
cance level 7 is equivalent to finding at least one null shape function in the 
set of all possible functions in the (1 — 7) confidence band. 

Proposition 1. Denote Qc = Wit) : a{t) = c, c G M} (the set of all real 
constant functions). The rejection rule of the hypothesis test for constant 
shape (Hq : a{t) € Qc) becomes 

ecr\CBj = where P{a{t) eCB^) = l--f. 

Using this rejection rule, the type I error is equal to 7. 

This proposition follows from the classical result on the equivalence of 
confidence intervals and hypothesis testing in Lehmann (1997). 

Space-time interaction. In our modeling approach, in order to account for 
the space-time interaction, we introduce an additional term in 7r(t, s) spec- 
ified by a set of effects i^r = {i^r,nm}n=i,...,N,m=i,...,M- The testing procedure 
for space-time interaction of the regression coefficient for the rth predictor 
reduces to 

Ho ■.au^=0 vs. Hq : a,^^ > 0. 

The null hypothesis implies that the association between the rth predictor 
and the response is separable in time and space. 

Although there are several approaches for testing the null hypothesis of 
zero variance component in linear mixed effects models, hypothesis tests 
that apply under multiple variance components have been investigated only 
recently. To test for space-time interaction under a multiple predictor model, 
we therefore use the approximations to the finite sample null distribution of 
the RLRT statistic in Greven et al. (2008). 
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3. Multilevel varying coefficient model. In this section we discuss the ex- 
tension of the varying-coefficient model in Section 2 to data with an intrinsic 
multilevel structure. 

3.1. The model. The observed data for the pth category (e.g., service 
provider) for p = 1, . . . , P are (Yp^ij, Xr^ij), where Yp^ij = Yp{ti,Sj) is a gener- 
alized response variable and X^^ij = Xr {ti , Sj ) the rth covariate observed at 
location Sj = {sji,Sj2) and time ti with E[l^^jj|X] = "fipiti, Sj)Xi^ij -|- • ■ • -|- 
lRp{tii Sj)Xji^ij where "frpit, s) for r = 1, . . . ,R are smooth coefficient func- 
tions. In our application, {i^,ij}j=i,...,Tj=i,...,5 stre the measures of service 
accessibility to the pth. service provider sites. 

To assess the association deviations of each of the kth group of processes 
from the global association pattern, we further decompose the regression 
coefficients as follows: 

(3.1) ^rpitii Sj) — ^ritii Sj) + T^rpi^i, Sj), 

where 7r(t, s) specifies the global association patterns and T]rp(ti,Sj) spec- 
ifies the group-specific deviations from the global association patterns. We 
further assume that the global effects 7r(t, s) take an additive form 

Jr(t,s) = Qr(t) +/3r(s), 

where the time- and space-varying regression coefficients follow the decom- 
position in (2.1). We also assume that the group-specific regression coeffi- 
cients are decomposed according to 

AI N 

Tjj-p(t, s) — Orp(t) -\- I3j.p(^s) + ^ ^ ^ ^ ^r,p,nmKtenipit ^p,'rn)-^sp(,S l^p^n)^ 

m=l n=l 

where arp{t) and f3rp{s) are decomposed using the radial spline basis simi- 
larly to the formulas in (2.1). We denote Kp^m, m = 1, . . . , M, the temporal 

knots used in the decomposition of the time-varying regression coefficient 

(s) 

and Kp^n, n = 1,. . . ,N, the spatial knots used in the decomposition of the 
space- varying regression coefficient for the pth service provider. For example, 
the decomposition of the regression coefficient arp{t) is 

M 

Orp(t) — Trpfi ~l~ Trp,lt ~l~ ^ ^ '^rp,m-^tcrapi\i ^p,ml)' 

m=l 

3.2. Estimation. Similar to the varying coefficient model in Section 2, 
we estimate the parameters in the multilevel varying coefficient model using 
the mixed effects model equivalence, resulting in a multilevel mixed effects 
model. 
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For the multilevel model, we need to impose a series of constraints on 
the fixed effects and on the selection of the temporal and spatial knots. For 
r = l,...,R, 

P P 

Trpfl = and y^^Trp,! = 0, 

p=l p=l 

p p p 

Srpfi = 0, ^ 5rp,l = and ^ 5rp,12 = 0. 

p=l p=l p=l 

Proposition 2. If the temporal and spatial knots are selected such that 

\k(t) _ JT) I ,(t) 

\'^mi,p '^m.2,p'i ^ " 
/or any mi,m2 € {1, ... , M}, and for any p,p' = 0,1, . . . , P {p^p'), 

ll'^ni,p ^n2,p'll ^ " 

for any ni,n2 G {1, . . . , A^}, and for any p,p' = 0,1, P {p ^ p'), where 
d(^) and d(^) are away from zero, then the model parameters in the multilevel 
model decomposition in Section 3.1 are identifiable. 

The proof of this proposition is provided in Supplemental Material 3 of 
this paper. 

3.3. Inference. Since making inference under the multilevel model pre- 
sented in the previous section implies making inference over all groups 
jointly, we need to correct for multiplicity. For instance, given that we need 
to evaluate the shape of the temporal global effects in the decomposition 
of the time- varying coefficients corresponding to the rth predictor, arp{t) 
for p = 1, . . . , P, we test multiple hypotheses [e.g., Hop-.arpit) constant] si- 
multaneously. For a small number of groups (P small), we can modify the 
approach discussed in Section 2.3 to account for the joint inference. Conse- 
quently, we estimate joint confidence bands: 

P{arp£CBrp,p = l,...,P) >l-p 

by correcting the confidence level of individual confidence bands for multiple 
inference using a Bonferroni correction; that is, estimate 1 — p/P confidence 
bands. Under the classical definition of the type I error for joint inference, 
we find that the test using 1 — p joint confidence bands is p, 

p 

type I error = ^ Prj/,.p(Gconst n CBrp = 0) < p. 
p=i 

Note that this correction will provide overly conservative confidence band 
estimates when P is large. 
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4. Case study: The equity of financial services. We proceed with the 
appHcation of the varying coefficient models to assess whether there are 
systematic disparities in the service accessibility with respect to various de- 
mographic variables. We focus on the equity of financial service accessibility 
in the state of Georgia over a period of 13 years, 1996-2008. We start with 
a description of the accessibility data followed by a brief exploratory analy- 
sis of the demographic variables. We continue with the presentation of the 
findings from the application of the varying coefficient models. 

4.1. Accessibility data. The site location data in this study were acquired 
from the Federal Deposit Insurance Corporation (FDIC). In our study we use 
data starting from 1996 to 2008. We geocoded the site location addresses 
using ArcGIS (ESRI) to obtain the service point locations in the service 
network: S = {si, . . . , s„} (n = 2,849 for Georgia). 

In service research the distance between a service site and its customers is 
commonly evaluated using the Euclidean or the Manhattan distance between 
the centroid of the neighborhood and the location of the closest service site. 
CIS road network data allows including more realistic route distances. For 
example, Talen (1997, 2001) uses the street-network distance to compute 
the distance between the centroid of the neighborhood and the site location. 
Lovett et al. (2002) use road distance and travel time by car. We acquired 
highway data for the whole U.S. (courtesy of the CIS Center at Georgia 
Institute of Technology) as well as a TIGER street-detailed network for 
Georgia and we took the average of the travel distances computed using 
both networks to obtain the distances d{ui,,Si;t) for b = l,...,B sample 
locations within a community and si,...,Sn service sites. Notably, none of 
the two networks provide highly accurate travel distances; therefore, the 
average over the distances computed using the two networks will provide 
more robust distance estimates. Finally, the travel cost C{ub,t) is computed 
as the average over the smallest three distances in {d{ub, Si;t),i = 1, . . . ,n}. 

Last, we obtain the utilization adjustment weights using the popula- 
tion counts acquired from the Environmental Systems Research Institute 
(ESRI). We use kernel smoothing [Diggle (1985)] to estimate the rate of 
point spatial processes. Using this approach, we obtain the population and 
service rate estimates at the sample locations, P{uh,t),b = 1,...,B and 
R{ub, t),b = 1, . . . , B, in year t. Further, we compute the utilization weights 
using W{uh,t) = P{ub,t)/R{uh,t),b = 1,. . . ,B, and along with the travel 
cost C{uh,t), we can finally obtain the accessibility measure at the com- 
munity level using equation (1.1). We apply this estimation procedure for 
all communities in Georgia and obtain the accessibility process Y{Us,t) for 
s = 1, . . . , S {S = 1,624) and t = 1996, . . . , 2008. In this research, census tracts 
are used as proxy for communities. According to the Census Bureau, census 
tracts are delineated with local input and intended to represent neighbor- 
hoods. 
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Fig. 2. (a) Logtravel cost — 1996; (b) Log travel cost — 2008; (c) Log accessibility — 1996; 
and (d) Log accessibility — 2008. 

Remark. Since the accessibility measure is an adjusted travel cost, we 
interpret it as follows. Large values of the travel cost or large values of the 
measure correspond to low accessibility to the service network. Therefore, 
if the measure values are, for example, increasing, the access to service is 
decreasing. Moreover, if the association of a demographic variable to the 
proposed measure is high, we infer that there is low accessibility with respect 
to the demographic variable. 

In the following discussion, we contrast "horizontal equity" [Figure 2(a), 
(b)] measured using the travel cost without adjusting for the "utilization" 
of a service operation [in equation (1.1), W{u,t) = 1 for any location u and 
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time point t] to ^^vertical equity" [Figure 2(c), (d)] which accounts for the 
expected utihzation of a service. 

Although difficult to assess visually, there are more extensive areas with 
lower (unadjusted) travel cost in 2008 than in 1996 in Georgia. On the other 
hand, the access to financial services is slightly lower in 2008 than in 1996 for 
highly populated regions, more specifically, Atlanta (see Supplemental Mate- 
rial 6 for the travel cost maps of metropolitan Atlanta). The primary reason 
for this contrast is that the increase in the number of new financial sites 
has a lower slope than the population growth in highly populated regions 
in Georgia. Consequently, these regions have weaker access compared to low 
density population areas, although the travel cost is small. These findings 
point to potential business opportunities for financial service providers. 

This comparison between travel cost without and with utilization-adjust- 
ment motivates the need for correcting the travel cost for the expected uti- 
lization of a service. In our subsequent analysis, we will only focus on the 
utilization- adjusted travel cost. 

4.2. Demographics data. In this study the demographic variables used 
to predict service accessibility include median household income, race and 
ethnicity data which are acquired from the Environmental Systems Research 
Institute (ESRI). The description of the methodology employed to obtain the 
demographic estimates at the census tract level is provided in Supplemental 
Material 5 of this paper. One has to bear in mind that the demographic 
estimates are measured with error which, in turn, will impact the estimates 
of the association between accessibility and the demographic variables. 

Since the boundaries of census tracts are updated by the Census Bureau 
every ten years, our data set includes a change of boundaries. The Census 
Bureau provides the so-called "relationship files" to document the revisions 
of the 1990 to 2000 census tract boundaries. We map the data collected 
before 1999 to 2000 boundaries using the information in these relationship 
files. 

Figure 3 shows the income level on the log scale and the percentages of 
Black, Hispanic and White populations for the last year of demographic data 
in this study. We do not show the plot of Asian percentages since overall in 
Georgia (except a small area in Atlanta) the percentage of Asian population 
is very low (close to zero). Contrasting the plots displaying the percentage of 
Black and White populations, we note that areas of high Black population 
have low White population and vice versa, pointing to significant segregation 
between black and white populations in Georgia. Indeed, the correlation is as 
high as —0.98, which suggests high collinearity between these two variables. 
At the same time, the collinearity between any other two demographic vari- 
ables is low (see Table 1). Since there is high collinearity only between White 
and Black populations, we fit the varying coefficient models separately for 
{income, Hispanic, White} and for {income, Hispanic, Black}. 
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Fig. 3. Demographic variables in Georgia. 



Table 1 

Correlation between demographic variables 



Income 


White 


Black 


Hispanic 


1.00 


0.18 


-0.27 


0.37 


0.18 


1.00 


-0.98 


-0.01 


-0.27 


-0.98 


1.00 


-0.16 


0.37 


-0.01 


-0.16 


1.00 



4.3. Varying coefficient model: Motivation, results and findings. In this 
paper we introduce a framework for studying the equity of service accessibil- 
ity for different groups of people identified by location or underlying socioe- 
conomic variables. The data consist of a series of maps characterizing the 
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access to financial services and a series of maps describing the demographic 
composition at the neighborhood level varying in time. The objective is to 
assess geographically-varying association patterns between accessibility and 
demographic variables over a period of several years. 

Simple visual inspection of a large number of maps (13 for the accessibility 
measure and 13 x 4 for four demographic variables) observed over a large 
geographic space goes beyond feasibility. Moreover, the existing models will 
only allow partial understanding of the dynamics in the equity of service 
accessibility. For example, space-varying coefficient models provide a one- 
year snapshot of the equity in service accessibility but will neither explain 
how it has changed over time nor account for the interaction between space 
and time associations. To model the space-time dynamics in the association 
between accessibility and socioeconomic variables jointly, we therefore apply 
the space-time varying coefficient model in Section 2. 

In the application of the space-time varying coefficient model to the data 
in this study, we selected a small number of temporal basis functions (M = 7) 
since we have a small number of time points; the space-time varying coef- 
ficients do not change significantly for various values of M. However, the 
estimated space-time varying coefficients vary with the number of spatial 
basis functions, A^. For small A^, the space-varying coefficients are smooth. 
Ruppert (2002) empirically suggests that after a minimum number of knots 
has been reached, the modeling bias is small. Therefore, we can control the 
modeling bias by using a large enough N] in our application N can be as 
large as 5 = 1,624. In contrast, the larger is, the more expensive the 
computation is. Consequently, we need to select N for an optimal trade-off 
between modeling bias and computational feasibility. To select N , we used 
a residual-based analysis suggested by Wood (2006). 

When interpreting the varying regression coefficients, one has to bear 
in mind that large values of the accessibility measure (population-adjusted 
travel cost) correspond to weak access to financial services. Moreover, sig- 
nificant association between accessibility and a demographic variable sug- 
gests that access to financial services is driven in part by the presence or 
the absence of the population group identified by the corresponding vari- 
able. 

In this section we summarize our findings based on Figures 4, 5, 6 and 7, 
which include the following: 

(1) The time-varying coefficients for income, % of Black population, % 
of Hispanic population and % of White population. 

(2) The space-varying association patterns for the four demographic 
covariates in 2008 calculated from 7(2008, s) = a(2008) + (3{s) + 
Interaction(2008,s). 

(3) The point locations of inequities with respect to the four demographic 
covariates in 2008. 
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(g) Accessibility vs. Hispanic in 2008 (h) Accessibility vs. White in 2008 

Fig. 4. Georgia; The time-varying coefficients (a) -(d) and the spatial relationship pat- 
tern m 2008 for four demographic covariates — income, percentage of Black, Hispanic and 
White populations. 

The output figures summarize the space-time relationships between ac- 
cessibility and the socioeconomic variables considered in this study. We high- 
light that without a rigorous modeling procedure, we cannot evaluate the 
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Fig. 5. Georgia; Inequity locations in 2008 with respect to (wrt) four demographic co- 
variates — income, percentage of Black, Hispanic and White populations. 

significance of tlie associations to service accessibility. Therefore, using the 
space-time varying model is important not only for estimation of these as- 
sociations but also for inference about their significance as described below. 

We define locations of inequity with respect to income to be the spatial 
units s such that 7incomc(2008, s) is statistically significantly positive (pos- 
itive correlation between income and utilization-adjusted travel cost). We 
also define locations of inequity with respect to race/ethnicity (percentage 
of Black, Hispanic and White populations) to be the spatial units s such 
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(g) Accessibility vs. Hispanic in 200S (h) Accessibility vs. White in 2008 

Fig. 6. Atlanta; The time-varying coefficients (a) -(d) and the spatial relationship pat- 
tern in 2008 for four demographic covariates — income, percentage of Black, Hispanic and 
White populations. 

that 7cthnicity(2008, s) is statistically significantly negative (negative corre- 
lation between percentage and utilization-adjusted travel cost). Statistical 
significance of the coefficients is derived from the simultaneous confidence 
bands of the spatial coefficients. Specifically, the coefficient at location s 
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Fig. 7. Atlanta; Inequity locations m 2008 with respect to (wrt) four demographic co- 
variates — income, percentage of Black, Hispanic and White populations. 

is statistically significantly positive with 95% significance level if the lower 
bound of the confidence interval at s is positive and it is statistically signifi- 
cantly negative if the upper bound of the confidence interval at s is negative. 

The time- varying coefficients corresponding to each demographic covari- 
ate and their confidence bands are in Figure 4(a)-(d) and Figure 6(a)-(d). 
Using the approach for evaluating the shape of the varying coefficients in 
Section 2.3, we infer that at the significance level of 95%, the time- varying 
coefficients for income are nonlinear, for the % of Black and Hispanic pop- 
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Table 2 

P-values for testing the significance of the space-time interaction in the varying 
coefficients for four demographic variables. 



Covariate 


Income 


% of Black 


% of Hispanic 


% of White 


p- value 


0.0018 


0.0009 


f=sO 


f»0 



ulations are linear whereas for the % of White population are constant. All 
space- varying coefficients are nonlinear. 

We also evaluate the significance of the interaction terms for all four 
demographic variables. We apply the testing procedure for the space-time 
interaction term described in Section 2.3. The p-values are provided in Ta- 
ble 2. For all four demographic variables, the space-time interaction terms 
are highly significant, which implies that the space-time interactions in the 
association patterns are statistically significant; this suggests that the inter- 
action term has a significant contribution to the spatial association patterns. 
Therefore, using a varying coefficient model with a space-time interaction 
term significantly contributes to a more accurate association analysis. 

Time-varying association patterns. Following the inference procedure for 
shape evaluation, we infer that over the past 13 years in the state of Georgia, 
the association between the access to financial services and two demographic 
variables, income and the percentage of Hispanic population, has strength- 
ened over time with a brief decrease in the last years, whereas the association 
between the access to financial services and the percentage of Black pop- 
ulation has weakened. This suggests that access to financial services has 
become more and more dependent on the income level of the residents in 
a community and whether they are of Hispanic descent but less dependent 
on the race of the population. 

Importantly, we cannot make inference about the magnitude of the associ- 
ation patterns since the constants for the temporal and spatial coefficients, tq 
and 5o, are nonidentifiable. For inference on the level and the direction (nega- 
tive or positive) of the service accessibility association we need to investigate 
the space-varying association year by year as discussed below. 

Space-varying association patterns. Controlling for race and ethnicity, 
the association between access to financial services and income level varies 
throughout the state of Georgia, with primarily weak positive association in 
the north but negative association in the south [Figure 4(e)]. This pattern 
is consistent with the income map in Figure 3(a); the income is consis- 
tently low in south and middle Georgia except for a few urban areas. This 
suggests that regions with low income population tend to also have lower 
access to financial services regardless of race and ethnicity. Moreover, there 
are only a few locations with statistically significant positive association be- 
tween income and utilization- adjusted travel cost [Figure 5(a)]. This implies 
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that although south Georgia consists primarily of low income population 
whereas north Georgia is more mixed with higher income population than 
south, financial services are present in both. 

The map of the Black population percentage is not as uniform as its 
association to accessibility; there is a high density of the Black population in 
south Atlanta and in mid to south Georgia but not in the north [Figure 3(b)] . 
On the other hand, there are several locations with statistically significant 
inequities as shown in Figure 5(b), although the association of the % of 
Black population to financial service access is weak in Georgia except for 
the upper north. Most of these locations are in urban areas. We therefore 
conclude that the inequities in access to financial services with respect to 
the Black population are present but low throughout Georgia. 

The association between travel cost and the % of Hispanic population is 
neither uniformly positive nor high [Figure 4(g)]. Areas of high Hispanic den- 
sity population have low but statistically significantly positive association 
[Figure 5(c)]. This indicates that the presence of financial services decreases 
with the increase in Hispanic population. 

In contrast, the association for White population is consistently weak 
throughout Georgia, although the White population density is high in most 
of Georgia except in the middle [Figure 3(d)]. Moreover, there are much 
fewer inequity locations than for Black and Hispanic populations and most 
are in rural areas [Figure 5(b)~(d)]. 

Since Atlanta is the largest city in Georgia with mixed income population 
and with a high percentage of Black, Hispanic and White populations, we 
applied the modeling procedures proposed in this paper to evaluate potential 
inequities in the Atlanta area and its surroundings. 

The only time-varying coefficient that changes its shape in comparison 
to Georgia is for the White population; it has a nonlinear shape. There is 
an increase in the impact of the % of White population on the access to 
services (equivalently, a decrease in the impact on the travel cost) up to 
2000 followed by a slower decrease thereafter. 

Significant inequities in the Atlanta area are with respect to income and 
the Black population. The association of the % of Black population to ser- 
vice access is negative and strong in many communities in south and north 
Atlanta, implying significant inequities even after controlling for the income 
level. However, the association map does not fully overlap with the density 
of the Black population; that is, while South Atlanta has a large Black pop- 
ulation [Figure 4(b), Supplemental Material 6], most inequities are in North 
Atlanta. 

There is a positive association between income and access to financial 
services in south Atlanta, an area with a predominantly low income pop- 
ulation [Figure 4(a), Supplemental Material 6]. Moreover, there is a weak 
association in north Atlanta and negative association in the east and west 



A SPACE-TIME VARYING COEFFICIENT MODEL 



23 



borders (possibly over-served areas). Therefore, when comparing the associ- 
ation map and its statistical significance to the map of the per capita income 
[Figure 4(a), Supplemental Material 6], we conclude that many communities 
with low and median income in central Atlanta have low access to financial 
services, as there is a statistically significant association between utilization- 
adjusted travel cost and income in these communities. 

The association pattern for the % of Hispanic and % of White population 
is largely negative, with just a few communities with statistically significant 
positive association, although both population groups are well represented in 
Atlanta [Figure 4(c), (d). Supplemental Material 6], indicating insignificant 
inequities with respect to the Hispanic and White populations in Atlanta. 

4.4. Multilevel varying coefficient model: Motivation, results and findings. 
Since monopoly is not common in service distribution, there is an intrinsic 
multilevel structure to service accessibility. At a higher level, we estimate 
the association to service accessibility over all service providers, whereas 
at a lower level, we estimate the deviations from the overall patterns as- 
sociated with each service provider in the network. One simple approach 
would be to apply the space-time varying coefficient model to the acces- 
sibility measure computed for each service provider separately. However, 
this approach only takes into account the variability within the network of 
each service provider but not the variability between service providers; in 
other words, this simple approach does not allow estimating the deviations 
from the between-providers association patterns. The multilevel varying co- 
efficient model in Section 3 estimates both overall space-time association 
patterns and the deviations from the overall pattern corresponding to each 
service provider. 

In this section we discuss the association accessibility patterns to five 
financial service providers: Bank of America (BoA), Branch Banking and 
Trust Company (BB&T), Regions Bank, SunTrust Bank and Wachovia 
(Wells Fargo since 2008). In Supplemental Material 7 of this paper, we in- 
clude the corresponding association patterns derived from the application 
of the multilevel space-time varying coefficient models. 

All five banks are in the top 10 largest banks in the U.S., with a variety 
of financial services including retail and commercial banking, mortgages, 
insurance products, trust services and securities brokerage. SunTrust and 
Regions banks are mainly based in southern states, BB&T is a national 
bank, whereas BoA and Wachovia are national banks with international 
subsidiaries. Although Bank of America has dominated the financial service 
market for many years, due to mergers and acquisitions, SunTrust Bank 
and Wachovia (Wells Fargo since 2008) Bank now dominate the market. 
The only bank that has not increased the number of brunches in Georgia 
and, in fact, has closed some of them, is Bank of America. 
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Time-varying association patterns. Following the inference procedure for 
shape evaluation, we infer that the time-varying deviations from overall as- 
sociation patterns denoted in this paper by r]rp{ti,Sj), where r is the index 
for the service provider and p is the predictor index, are all approximately 
zero (not statistically significant) except for the deviation coefficients corre- 
sponding to income. Therefore, over the past 13 years in the state of Georgia, 
the association between the access to financial services and income is posi- 
tive and has strengthened for BB&T and Bank of America beyond the global 
upward trend, indicating stronger association between service access and the 
income level in 2008 as compared to 1996. The income deviation coefficient 
for Wachovia is decreasing over time but negative, implying a decrease in 
association with respect to income. Finally, there are not significant system- 
atic disparities in the provider-specific deviations from the overall association 
patterns for the ethnicity and race demographic variables. 

Space-varying association patterns. The lowest association between ser- 
vice access and the four demographic variables is for Bank of America, Re- 
gions and Wachovia. The accessibility association patterns for these three 
banks do not deviate significantly from the global trends. The association 
of the percentage of Black and White populations to service accessibility is 
approximately zero, whereas the association to the percentage of Hispanic 
population is weak with mixed association throughout Georgia. 

In urban and rural Georgia, the association between the access to BB&T 
services and the income level is highly positive, whereas the association for 
the % of Black and Hispanic populations is consistently weak and for the % 
of White population is approximately zero. Notably, the most significant 
inequities in 2008 for BB&T are with respect to income and they have 
increased over the 13 year period. For SunTrust, the service accessibility 
association to income is strong in south Georgia. There is weaker association 
in urban areas than in rural areas. Moreover, the association of the % of 
Black and Hispanic populations is weak. 

The most significant inequities in Atlanta are with respect to income — 
high association between income and service accessibility for Bank of Amer- 
ica and Regions in south Atlanta, for BB&T and SunTrust throughout At- 
lanta. There are contrasting association patterns for north and south At- 
lanta which also differ in their demographic decomposition. Generally, the 
accessibility association is low for the White population after controlling 
for income; this suggests equitable accessibility to financial services for the 
White population. 

5. Final considerations. The methodological contributions described in 
this paper are twofold. First, we introduce a framework for the study of the 
equity of service accessibility across population groups with various demo- 
graphic characteristics. This study allows characterization of the geograph- 
ically varying equity patterns over a period of several years. Second, we 
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investigate spatio-temporal estimation methods, which use the underlying 
structure of varying coefficient models. The first model estimates space-time 
varying association to a response variable (e.g., accessibility measure) of a se- 
ries of predictors (e.g., demographic variables) jointly. The second model 
extends the first model to a response variable with a multilevel structure. 
Because of the complexity of the model parameters, we propose a simplified 
inference procedure based on confidence bands which allows evaluation of 
the shape of the varying coefficients. 

We note that different service accessibility measures will provide different 
accessibility maps, and, therefore, different conclusions will be drawn for the 
study of service accessibility equity. In this paper, the underlying measure 
is defined as the utilization- adjusted travel cost; in Section 4.1 we compared 
the accessibility maps with and without correction for utilization and we 
concluded that the two measures will provide different perspectives in the 
equity of service accessibility, primarily for areas with low density popula- 
tion. A comparison study of accessibility measures is beyond the scope of 
this paper. 

From the analysis of service accessibility using the space-time varying 
coefficient model, we identified significant but low inequities in some regions 
of Georgia with respect to income after controlling for race and ethnicity, 
and with respect to Black and Hispanic populations after controlling for 
income. These inequities have increased over time. The most predominant 
inequities in Atlanta are for Black population, although they have decreased 
over time. The association between income and service access is largely pos- 
itive in Atlanta, suggesting potential inequities with respect to income as 
well. After accounting for service utilization, there are more significant in- 
equities in urban areas than in rural areas; this may be due to the fact that 
the population in rural areas is more homogeneous. 

In the analysis of service accessibility using the multilevel space-time 
model, we found for Georgia, and particularly for Atlanta, that the devi- 
ations from the between-provider association patterns are very insightful. 
Specifically, we learn, for example, that income-driven inequities for BB&T 
are significantly stronger when contrasted to the overall association to acces- 
sibility, whereas for other service providers, for example, Bank of America 
and Wachovia (Wells Fargo), there are not significant deviations from the 
between-providers association patterns. 

Importantly, one challenge in space-time varying coefficient model esti- 
mation is whether an assumed pattern in a multiple predictor model can 
actually be recovered. For this, we conducted a simulation study with two 
predictors. The association patterns for both predictors are nonseparable in 
space and time, a more realistic simulation framework. For this simulation, 
the estimated coefficients are accurately estimated. We also evaluated the 
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coverage and the power of the shape evaluation procedure discussed in Sec- 
tion 2.3. The power is lower for the time-varying components than for the 
space-varying ones; in this simulation study, the number of spatial points is 
S = 300, whereas the number of time points is T = 15, which may lead to 
lower accuracy in the shape evaluation of the varying coefficients. Moreover, 
the power depends on how smooth and close to the null hypothesis the shape 
function is. 

We note that different service accessibility measures will provide different 
accessibility maps, and, therefore, different conclusions will be drawn for the 
study of service accessibility equity. In this paper the underlying measure is 
defined as the utilization- adjusted travel cost; in Section 4.1 we compared 
the accessibility maps with and without correction for utilization and we 
concluded that the two measures will provide different perspectives in the 
equity of service accessibility. A comparison study of various accessibility 
measures is beyond the scope of this paper. 

One limitation of the study of service accessibility equity for large ge- 
ographic regions (e.g., the US) using the space- varying coefficient models 
introduced in this paper is the computational aspect. Although we have re- 
duced the estimation of multiple space-time varying coefficients with differ- 
ent smoothing levels to a simple mixed effects model, the estimation prob- 
lem remains computationally intensive. A large geographic space requires 
a large number of knots, which in turn results in a model with a large 
number of random effects. In addition, because of the model decomposition 
into separable and nonseparable space-time coefficients, the number of vari- 
ance components is large even for a small number of predictors (e.g., for 
three predictors in our study, we have a total of 12 variance components 
for the simple varying coefficient model but as many as 42 for the multi- 
level model). Ongoing research focuses on overcoming these challenges by 
using a backfitting estimation algorithm in the presence of multiple predic- 
tors. 
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Supplement 1. Varying-coefficient model — Decomposition of the design 
matrix under the tensor-product decomposition of the space-time varying 
coefficients. 

Supplement 2. Varying-coefficient model — Derivation of the confidence 
bands for the space and time varying coefficients. 

Supplement 3. Varying-coefficient model — A simulation study under mul- 
tiple predictors. 

Supplement 4. Varying-coefficient model — Proof of Proposition 2. 
Supplement 5. Case study — Description of ESRI data. 
Supplement 6. Case study — Accessibility maps for Atlanta area. 
Supplement 7. Case study — Results and maps for the provider-level ac- 
cessibility analysis. 
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